In this notebook we use QGIS to create a shareable terrain model with a data overlay, which can be shared on a web server, without typing a single line of code. Let's go!
To inetract with the maps we make below, check this notebook out in the ipython notebook viewer:
In this approach we download all our data from NCI to our machine (or could use it in the VDI), because we need to modify bands - which QGIS is not happy to do for web services (WMS or WCS).
For topography we need a terrain model in a raster format, e.g. GeoTIFF, which covers Canberra. We can head to the NCI Elevation collection here:
http://dapds00.nci.org.au/thredds/catalog/rr1/Elevation/catalog.html
...and look for SRTM 1 second elevation - good enough for this job. If you are happy with ESRI grids, navigate to the tile collection here:
http://dapds00.nci.org.au/thredds/catalog/rr1/Elevation/1secSRTM_DEMs_v1.0/DEM/Tiles/catalog.html
Data are organised in folders named by longitide and latitude of the south west (bottom left) corner. For Canberra we need 149, -36.
...and for now, this is probably the best method. We could make an OpenDAP or WCS request for a subset, but that would be coding! The pull of the dark side is strong - so here is a link that gets a GeoTIFF from the SRTM tile, using WCS:
Now, import the resulting GeoTIFF into QGIS as a raster data source:
Now zoom into the DEM so that it covers the entire display window, then head to the 'web' menu. Choose 'qgis2threejs'. In the resulting dialog box, click 'world' in the left pane, and find 'vertical exaggeration'. Set it to 10.
In 'DEM', the one active layer (your GeoTIFF) should be preselected. Click 'run', and you should get a web browser opening with a 3D model inside!
The expected output is shown below - try viewing the notebook here:
...you should be able to move the map, zoom in and out, and generally inspect a terrain model.
In [2]:
###ignore this block of code - it is required only to show the map in iPython - you won't need it!
from IPython.core.display import display, HTML
display(HTML('<iframe width="800" height="600" frameborder="1" scrolling ="no" src="./qgis2threejs/ACT_elevs_test_1.html"></iframe>'))
We might be interested in land cover attributes for our region, so let's get them! How about photosynthetic vegetation for 2015?
Using WCS again - click on the WCS link, and look for a
...so we put a question mark after the dataset name, then add the rest of the labels describing the thing we want afterward, in any order, separated by ampersands:
(woah. There was a glitch in the matrix - if we didn't write out that URL, you would have had to code just now).
Add the resulting GeoTIFF to QGIS as a raster data source, just like the DEM! Once you have a style you're happy with, soom to the desired extent and use the qgis2threejs plugin to make a new map.
In [3]:
display(HTML('<iframe width="800" height="600" frameborder="1" scrolling ="no" src="./qgis2threejs/act_elevs_plus_greenveg.html"></iframe>'))
Now we have an elevation map coloured by green vegetation! But that's only a pretty picture.
We will try to show:
...and no code, Only button clicking.
Back to QGIS. We need another plugin - zonal statistics! We know how to install plugins, so get to it. When ready, we'll make a hilliness proxy first.
Next we will make a vegetation cover proxy
..and finally, can we come up with a metric of vegetation cover as a function of block hilliness? A question here might be 'do land owners tend to clear flatter blocks more than hilly blocks?'. Can we answer it using web services data, QGIS, some clicks - and then make a very pretty, and interactive map?
On to the first question - how do we get ACT block data? Head to ACTMAPi and find ACT blocks. Here's a shortcut:
http://actmapi.actgov.opendata.arcgis.com/datasets/afa1d909a0ae427cb9c1963e0d2e80ca_4
Find the 'API' menu box, click the down arrow and then copy the URL inside the JSON option:
http://actmapi.actgov.opendata.arcgis.com/datasets/afa1d909a0ae427cb9c1963e0d2e80ca_4.geojson
In QGIS, add a new vector layer. Choose 'service' from the options, and paste the URL into the appropriate box:
Wait a while! Now because QGIS you'll need to save the result as a shapefile before we can do anything with it. Save the layer, close the GeoJSON layer and reopen your new ACT blocks shapefile.
...it's probably simpler to just download and open the shapefile - but now we know how to add vector layers from a web service> either way, the result is something like this:
We don't need to make the blocks pretty yet, We'll do something!
Using the SRTM DEM, let's make a proxy of block hilliness. Check that you have the QGIS zonal statistics plugin (Raster -> zonal statistics). If not, install it the same way you installed qgis2three.js. Once you have it, open the plugin. Choose the DEM as the raster layer, and use band 1. Then choose your ACT blocks vector layer. In the statistics to calculate, pick an appropriate set - but include standard deviation - this is our roughness proxy.
Run the plugin, then open the properties box of the ACT blocks layer and colour your blocks by standard deviation. The ZOnal statistics plugin has looped through all the polygons in the blocks layer and computed descriptive statistics of the underlying DEM. And there you have it - ACT blocks coloured by the standard deviation of the elevation they contain, as a proxy for hilliness.
Now lets make a cool map! zoom in so that you have a region you're happy with occupying your map view, and open the qgis2threejs plugin.
In the left pane, under polygon option choose your ACT blocks layer. Options for styling it appear on the right - which I'll recreate again because I crashed QGIS again at a bad time (Save, fella, save!)
What we see here are cadastral blocks, with our 'hilliness proxy' displayed by colour and extruded column height. Darker blue, and taller columns are hillier! Click on the map to inspect features, click and drag to move, scroll to zoom.
In [4]:
display(HTML('<iframe width="800" height="600" frameborder="1" scrolling ="no" src="./qgis2threejs/act_block_hilliness_proxy.html"></iframe>'))
We have some elevation data, we have some cadastral data, we have some data about photosynthetic vegetation cover. We can do our quick visualisation of whether block hilliness (as determined by SRTM height standard deviation for each block) is related to photosynthetic vegetation cover (as determined by the median of vegetation cover inside each block). We can set this up as follows:
In this scheme, if our hypothesis is that hillier blocks are less cleared, dark green blocks will visualise as taller columns. Lets test it out!
In [5]:
display(HTML('<iframe width="800" height="600" frameborder="1" scrolling ="no" src="./qgis2threejs/veg_mean_colours.html"></iframe>'))
In [ ]:
In [ ]: